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The  reliability  modeling  of  a  module  in  a  turbine  engine  requires  knowledge  of  its  failure  rate,  which  can  be  estimated  by  identifying 
statistical  distributions  describing  the  percentage  of  failure  per  component  within  the  turbine  module.  The  correct  definition  of  the 
failure  statistical  behavior  per  component  is  highly  dependent  on  the  engineer  skills  and  may  present  significant  discrepancies 
with  respect  to  the  historical  data.  There  is  no  formal  methodology  to  approach  this  problem  and  a  large  number  of  labor  hours 
are  spent  trying  to  reduce  the  discrepancy  by  manually  adjusting  the  distributions  parameters.  This  paper  addresses  this  problem 
and  provides  a  simulation-based  optimization  method  for  the  minimization  of  the  discrepancy  between  the  simulated  and  the 
historical  percentage  of  failures  for  turbine  engine  components.  The  proposed  methodology  optimizes  the  parameter  values  of 
the  component’s  failure  statistical  distributions  within  the  components  likelihood  confidence  bounds.  A  complete  testing  of  the 
proposed  method  is  performed  on  a  turbine  engine  case  study.  The  method  can  be  considered  as  a  decision-making  tool  for 
maintenance,  repair,  and  overhaul  companies  and  will  potentially  reduce  the  cost  of  labor  associated  to  finding  the  appropriate 
value  of  the  distribution  parameters  for  each  component/failure  mode  in  the  model  and  increase  the  accuracy  in  the  prediction  of 
the  mean  time  to  failures  (MTTF). 


1.  Introduction 

There  are  several  failure  modes  in  gas  turbine  engines.  Typ¬ 
ically,  the  failure  rates  of  different  components  are  recorded 
throughout  the  life  of  the  engine.  Upon  obtaining  the  failure 
rates  of  such  components,  the  first  step  is  to  conduct  a 
statistical  distribution  fitting  to  the  initial  data  per  component 
or  failure  mode.  The  mean  percentage  of  failures  computed 
while  considering  the  fitted  statistical  distribution  may  not 
be  similar  to  the  historical  data  due  to  the  lack  of  quantity  or 
quality  of  the  data.  Therefore,  a  common  practice  is  to  per¬ 
form  some  adjustments  in  the  distribution  parameters  based 
on  human  intelligence  and  the  experience  of  the  engineers. 
Currently,  there  is  no  formal  methodology  to  approach  this 
complex  problem. 

The  reliability  modeling  of  turbine  engines  is  a  complex 
stochastic  system.  The  complexities  that  arise  when  random¬ 
ness  is  embedded  within  a  system  make  the  analysis  of  these 


systems  a  difficult  task.  Unfortunately,  epistemic  uncertainty 
is  a  common  and  unavoidable  characteristic  among  real- 
world  systems.  The  development  of  simulation  as  means  to 
evaluate  stochastic  systems  enhances  the  ability  to  obtain  per¬ 
formance  measure  estimates  under  several  given  conditions. 
Furthermore,  these  performance  measures  are  much  more 
accurate  compared  to  the  estimations  of  analytical  tech¬ 
niques,  which  often  make  crude  assumptions  about  the 
systems  condition  and  operation. 

Computer  simulations  are  highly  effective  in  answering 
evaluative  questions  concerning  a  stochastic  system.  How¬ 
ever,  it  is  often  necessary  to  determine  the  values  of  the 
models  decision  parameters  such  that  a  performance  mea¬ 
sure  is  maximized  or  minimized.  This  type  of  problems  can 
be  tackled  by  simulation-based  optimization  methods  [1]. 

This  paper  aims  to  develop  a  simulation-based  optimiza¬ 
tion  method  to  reduce  the  discrepancy  between  simulated 
and  historical  failures,  which  will  not  only  significantly 
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Figure  1:  Series  system  for  the  mock  engine  model. 


reduce  the  cost  of  labor  hours  associated  to  performing  this 
activity  manually  but  will  also  find  an  appropriate  value  of  the 
distribution  parameters  for  each  turbine  component/failure 
mode  so  that  the  simulations  outcome  better  agrees  with  the 
real-life  data. 

This  research  was  done  in  collaboration  with  an  aerospace 
industry’s  maintenance,  repair,  and  overhaul  company 
referred  to  as  ABC  Company.  A  54-component  turbine 
engine  is  considered  as  a  case  study.  For  confidentiality 
reasons,  the  names  of  the  components,  distributions,  and 
distributions’  parameters  are  not  mentioned. 

The  structure  of  the  paper  is  organized  as  follows. 
Section  2  presents  the  proposed  model  and  its  notation. 
Section  3  discusses  the  details  of  the  proposed  methodology. 
In  Section  4,  the  case  study  is  presented  along  with  a  design  of 
experiments  (DOE)  to  tune  the  parameters  of  the  optimiza¬ 
tion  algorithm.  Section  5  shows  the  results  of  the  case  study. 
Finally,  concluding  remarks  and  future  work  are  presented  in 
Section  6. 

2.  Turbine  Engine  Model 

The  engine  model  is  developed  using  the  following  assump¬ 
tions:  (1)  the  component  events  are  independent  of  each 
other,  (2)  the  component  failure  rates  interact  in  a  series 
fashion  (see  Figure  1),  and  (3)  some  of  the  components 
may  have  accumulative  time  (e.g.,  an  item  that  has  been 
repaired  will  recover  its  previous  cumulative  time).  The 
model  considers  the  lowest  TTF  associated  to  any  component 
as  the  system’s  time  to  failure  (i.e.,  whichever  component  will 
have  the  lowest  TTF  at  each  run  of  the  simulation  will  cause 
the  system  to  fail) 

9  ( Tk )  =  min  (f1;  t2,  t3,  tn) .  (1) 

Equation  (1)  computes  the  system’s  TTF  where  g(Tk)  is 
the  system’s  TTF  and  t{  is  the  TTF  of  ith  component. 

3.  Methodology 

The  proposed  simulation-based  optimization  method 
searches  for  a  high  quality  solution  (i.e.,  the  value  of  the 
distribution’s  parameters  of  the  components  within  their 
feasible  regions)  with  the  use  of  the  simulated  annealing 
(SA)  metaheuristic.  The  proposed  method  is  broken  down 
into  two  phases. 

Phase  I  consists  of  a  Monte  Carlo  simulation  to  obtain 
the  simulated  percentage  of  failure  per  component,  given 
an  initial  set  of  distribution  parameters.  The  simulated 
percentages  of  failures  from  Phase  I  are  used  to  compute  the 
initial  discrepancy  between  the  simulated  and  the  historical 


percentages  of  failures,  which  is  expressed  as  the  sum  of 
squares  of  errors  (SSE).  Nevertheless,  other  combinations  of 
feasible  distribution’s  parameters  values  may  lead  to  lower 
SSE  values,  (i.e.,  other  high  quality  solutions  could  exist  in 
nearby  neighborhoods).  These  neighborhoods  are  explored 
in  the  next  phase  of  the  method. 

In  Phase  II,  a  SA-based  search  procedure  is  used  to 
explore  the  neighborhood  within  the  feasible  region  of  the 
distributions’  parameters  in  order  to  locate  solutions  of 
potentially  higher  quality.  The  objective  is  to  minimize  the 
SSE  obtained  from  Phase  I.  The  SA-based  procedure  works 
well  for  Phase  II  because  it  has  the  ability  to  move  away  from 
local  minima  in  contrast  with  gradient  methods  that  cannot 
overcome  local  minima.  The  SA  neighborhood  function  is 
defined  based  on  the  bounds  of  the  distribution  parameters  of 
each  component.  The  following  section  discusses  the  imple¬ 
mentation  in  further  detail. 

3.1.  Monte  Carlo  Simulation  to  Obtain  the  Failure  Time  of 
Each  Component.  A  Monte  Carlo  simulation  is  often  used  to 
obtain  not  only  the  simulated  failure  times  per  component, 
but  also  the  distribution  of  the  statistical  estimates  of  the 
system’s  TTFs:  MTTF  and  STTF. 

At  this  stage  of  the  proposed  methodology,  the  distribu¬ 
tion  and  the  distributions’  parameters  values  per  component 
(obtained  from  initial  data)  are  considered  as  primary  inputs. 
The  data  set  for  each  component  contains  suspensions  (cen¬ 
sored  data)  which  represent  the  units  that  have  not  failed 
by  the  failure  mode  in  question.  They  may  have  failed  by 
different  failure  modes  or  not  failed  at  all  [2]. 

In  the  case  of  life  data,  these  data  sets  are  composed  of 
units  that  did  not  fail.  For  example,  if  five  units  were  tested 
and  only  three  had  failed  by  the  end  of  the  test,  it  would  have 
right  censored  data  (or  suspension  data)  for  the  two  units  that 
did  not  failed.  The  term  right  censored  implies  that  the  event 
of  interest  (i.e.,  the  time  to  failure)  is  to  the  right  of  our  data 
point.  In  other  words,  if  the  units  were  to  keep  on  operating, 
the  failure  would  occur  to  the  right  on  the  time  scale  [3]. 

A  random  number,  drawn  from  a  uniform  distribution 
over  the  interval  [0, 1]  is  used  to  generate  a  failure  time  from 
the  distribution’s  inverse  cumulative  distribution  function 
(CDF).  For  each  component  of  the  model,  the  formula  for  the 
CDF  corresponding  to  the  distribution  of  each  component  is 
used  and  then  solved  for  t  (the  failure  time).  Table  1  displays 
the  equations  to  generate  failure  times  ( t ). 

For  the  normal  or  lognormal  distributions,  there  are  no 
closed-form  expressions  but  there  are  good  mathematical 
approximations  for  that  purpose.  The  approximation  utilized 
in  this  paper  is  the  one  developed  by  Weisstein  [4],  which 
considers  the  error  function  (erfc).  The  definition  for  the  erfc 
can  be  found  in  the  Appendix.  The  inverse  CDF  equations, 
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Table  1:  Inverse  CDF  equations  for  different  distributions  in  the  case  study. 


Distribution 

Inverse  CDF 

Normal 

-  V2<7erfcinv  j 

>(i-(<i -»(-(0-c  (-mm))} 

Fh 

i 

=L 

+ 

Log  normal 

e[-(V2Cerfcinv(2((-(l-F(0)(l-(0.5erfc(-((lnT-A)/V20))))+l)))+A]  _  j 

Weibull  (2  parameters) 

«((-ln((l-F(t))  <T(r/^))W)-T 

Exponential 

H-FO)] 

shown  in  Table  1,  include  cumulated  time  (T);  that  is,  as 
different  individual  components  are  replaced,  the  different 
operating  times  can  be  adjusted  for  the  components  whose 
reliability  has  been  restored  (through  either  refurbishment  or 
replacement).  Components  that  are  not  replaced  are  returned 
to  service  (i.e.,  their  operating  times  are  not  different  than 
when  they  were  removed  from  service).  This  reliability  model 
will  be  used  to  compute  the  systems  conditional  probability  of 
failure  once  it  has  been  repaired.  The  inverse  CDF  equations 
are  derived  from  the  conditional  reliability  equations.  If  there 
is  no  cumulated  time  (e.g.,  T  =  0),  then  the  inverse  CDF 
equations  would  be  transformed  into  the  unconditional  form. 

The  systems  TTF  is  computed  considering  the  compo¬ 
nent  with  the  lowest  TTF  as  the  system  failure  rate  (1).  Thus,  a 
single  iteration  of  the  loop  computes  a  single  system  TTF  for 
the  model.  Statistical  estimates  (mean  and  standard  deviation 
of  the  system  TTFs  MTTF  and  STTF,  resp.)  as  well  as  the 
simulated  percentages  of  failures  for  each  component  are 
obtained  after  performing  all  the  iterations  of  the  Monte 
Carlo  simulation. 

3.2.  Simulated  Annealing- Based  Optimization  Procedure.  The 
simulated  annealing  (SA)  algorithm  is  a  well-known  local 
search  metaheuristic  used  to  address  discrete,  continuous, 
and  multiobjective  optimization  problems.  Its  ease  of  imple¬ 
mentation,  convergence  properties,  and  ability  to  escape  local 
optima  have  made  it  a  popular  technique  during  the  past  two 
decades.  The  application  of  the  SA  algorithm  to  optimization 
problems  emerges  from  the  work  of  Kirkpatrick  et  al.  [5]  and 
Cerny  [6].  In  these  pioneering  works,  the  SA  was  applied  to 
graph  partitioning  and  very-large-scale  integration  design. 
In  the  1980s,  the  SA  algorithm  had  a  major  impact  on  the 
field  of  heuristic  search  for  its  simplicity  and  efficiency  in 
solving  combinatorial  optimization  problems.  Moreover,  the 
SA  algorithm  has  been  extended  to  deal  with  continuous 
optimization  problems  [7-9]. 

The  SA  algorithm  is  based  on  the  principles  of  statistical 
mechanics  whereby  the  annealing  process  requires  heating 
and  then  a  slow  cooling  of  a  substance  to  obtain  a  strong 
crystalline  structure.  The  strength  of  the  structure  depends 
on  the  rate  of  cooling  of  the  substance.  If  the  initial  tem¬ 
perature  is  not  sufficiently  high  or  a  fast  cooling  is  applied, 
imperfections  (metastable  states)  are  obtained.  In  this  case, 
the  cooling  solid  will  not  attain  thermal  equilibrium  at  each 
temperature.  Strong  crystals  are  grown  from  careful  and  slow 


cooling.  The  SA  algorithm  simulates  the  energy  changes  in 
a  system  subjected  to  a  cooling  process  until  it  converges  to 
an  equilibrium  state  (steady  frozen  state).  This  scheme  was 
developed  in  1953  by  Metropolis  et  al.  [10]. 

The  SA-based  procedure  has  been  successfully  imple¬ 
mented  to  address  problems  in  manufacturing  enterprises 
such  as  [11,  12].  The  interested  reader  is  referred  to  the 
works  of  Eglese  [13],  Aarts  and  Lenstra  [14]  and  Suman  and 
Kumar  [15],  for  a  comprehensive  description  of  the  simulated 
annealing  metaheuristic. 


3.2.1.  Optimization  Problem  Definition,  Objective,  and  Con¬ 
straints.  After  all  the  iterations  of  the  Monte  Carlo  simula¬ 
tion,  the  percentage  of  failures  is  obtained  for  each  compo¬ 
nent  and  it  is  compared  to  the  historical  data  to  obtain  the 
squared  error.  The  discrepancy  or  the  squared  error  between 
the  simulated  and  historical  percentage  of  failure  is  obtained 
for  each  component  to  conform  the  sum  of  squared  errors 
(SSE). 

The  objective  function  for  the  optimization  process  is 
shown  in  (2)  and  aims  to  minimize  the  sum  of  squared 
errors  (SSE)  of  the  systems  percentage  of  failures  subject  to 
constraints  which  are  the  confidence  bounds  (contour  plots) 
on  the  distributions  parameters: 


min 


n 

SSE  =  y  (Xi  -  xtf 

i=  1 


(2) 


where  x{  is  the  simulated  percentage  of  failures,  x{  is  the 
historical  percentage  of  failures,  and  n  is  the  number  of 
components. 

The  simulation-based  optimization  method  will  propose 
new  values  of  the  distribution  parameters  (P/;),  where  P/;- 
is  the  z'th  distribution  parameter  of  jth  component  within 
their  contour  plots  (see  Figure  2).  According  to  Wasserman 
[16],  the  likelihood  contour  plot  can  be  drawn  satisfying  the 
following  relationship: 

lnL{Pij)  =  lnL*  ~\x2l,a’  (3) 

where  In  L*  is  the  log-likelihood  of  the  jth  component  at  the 
estimated  distribution  parameters  using  MLE,  In  L(P/;)  is  the 
log-likelihood  of  the  jth  component  for  the  perturbed  zth 
distribution  parameter,  and  a -l-  confidence  level. 
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Figure  2:  Feasible  region  (contour  plot)  for  the  distribution  param¬ 
eters  of  one  component. 


In  case  of  unfeasibility,  that  is,  the  distributions  parame¬ 
ters  at  the  new  state  not  satisfying  (3),  one  of  the  distribution 
parameters  is  set  to  its  nearest  bound  (upper  or  lower  bound, 
whichever  is  closer)  and  solved  for  the  other  distribution 
parameter  so  that  they  satisfy  (3).  Equation  (4)  shows  the 
equation  of  likelihood  contour  plot  for  the  Weibull  two- 
parameter  distribution  as  an  example: 

In  L  (a  (lower  bound  or  upper  bound) ,  /3) 

1  0  W 

=  In  L* - y  1 

2  **  1  yCl 

The  shaded  region  in  Figure  2  shows  an  example  of 
the  parameters  feasible  region  for  the  simulation-based 
optimization  method. 

A  schematic  diagram  of  the  SA-based  procedure,  used  to 
determine  the  best  value  of  the  distributions  parameters,  is 
presented  in  Figure  3. 

Obtaining  the  confidence  bounds  of  the  distributions 
parameters  is  an  integral  part  of  the  proposed  methodology 
There  are  multiple  methods  for  obtaining  the  bounds  of  these 
parameters.  The  following  section  provides  the  necessary 
background  behind  the  selection  of  the  likelihood  ratio  based 
on  the  confidence  bounds  as  the  preferred  method  as  well 
as  the  steps  for  obtaining  the  bounds  of  the  distributions 
parameters. 

3.2.2.  Confidence  Bounds  on  the  Distribution  Parameters. 
Our  data  set  contains  suspensions  (censored  data);  therefore, 
the  procedure  for  computing  confidence  bounds  includes 
censored  data  in  the  parameter  estimation.  In  this  section,  we 
present  the  methodology  to  obtain  the  likelihood  ratio-based 
confidence  bounds  of  the  distributions  parameters. 

In  the  general  case,  the  distribution  parameters  are 
estimated  using  the  maximum  likelihood  estimation  (MLE) 
method  with  a  modified  likelihood  function  for  censored 
data.  Using  these  estimated  parameters,  confidence  bounds 
can  be  calculated.  The  parameters  are  estimated  using 


“mlecustom”  function  of  MATLAB.  The  general  form  of  the 
likelihood  function  considering  censored  samples  is  [2] 

L  =  fl/(^)fl(1-F(r;))>  W 

i= 1  j= 1 

where  r  is  the  number  of  units  run  to  failure,  k  is  the  number 
of  unfailed  units,  xl9  x2,  x3, . . . ,  xr  are  the  known  failure 
times,  and  Tl9T2>T39. . .  yTk  are  the  operating  times  on  each 
unfailed  unit. 

For  example,  if  the  time  to  failure  distribution  is  Weibull- 
distributed,  the  modified  likelihood  function  for  censored 
samples  would  be 


The  “mlecustom”  function  of  MATLAB  was  used  to 
compute  the  bounds  of  the  distributions  parameters  using  the 
normal  approximation  method.  The  normal  approximation 
method  for  confidence  interval  estimation  is  used  in  most 
commercial  statistical  packages  because  of  the  relative  easi¬ 
ness  for  the  bounds  computation.  However,  the  performance 
of  such  procedures  could  be  poor  when  the  sample  size  is  not 
large,  or  when  heavy  censoring  is  considered  [17]. 

On  one  hand,  since  the  input  data  set  contains  heavy  cen¬ 
soring,  the  computation  of  confidence  bounds  on  parameters 
using  the  normal  approximation  is  not  recommended.  On  the 
other  hand,  Fisher  matrix  bounds  tend  to  be  more  optimistic 
than  the  nonparametric  rank  based  bounds.  This  may  be  a 
concern,  particularly  when  dealing  with  small  sample  sizes. 
Fisher  matrix  bounds  are  too  optimistic  when  dealing  with 
small  sample  sizes  and,  usually,  it  is  preferred  to  use  other 
techniques  for  calculating  confidence  bounds,  such  as  the 
likelihood  ratio  bounds  [2] .  Hence,  the  likelihood  ratio-based 
confidence  bounds  estimation  method  is  preferred. 

Without  loss  of  generality,  the  use  of  the  likelihood 
ratio  statistic  for  developing  confidence  intervals  about  a 
parameter  of  interest  (y)  can  be  described.  The  likelihood 
ratio  (LR)  test  statistic  provides  a  formal  framework  for 
testing  the  hypothesis:  H0:  y  =  y0  and  Hx:y  ±  y0. 

As  its  name  implies,  the  LR  statistic  test  is  a  ratio  of 
likelihood  functions.  However,  it  is  more  convenient  to  work 
with  the  log  form,  which  is  computed  as  the  difference 
between  two  log-likelihood  expressions.  Specifically,  H0  is 
rejected  at  some  level  of  confidence  (1  -  a),  if 

2(lnL* -InL*  (y0))  >  yifl)  (7) 

where  InL*  is  the  log-likelihood  function  evaluated  at  the 
maximum  likelihood  (ML)  point  In  L*  (/3,  y)  and  In  L*  (y0)  = 
max^ln  L(y0,/3).  If  /?(y0)  is  denoted  by  the  solution  to 
max^  In  L(y0,  /3),  then 

lnL*(y0)  =  lnL(y0,^(y0)).  (8) 

An  asymptotically  correct  confidence  interval  or  region 
on  y  consists  of  all  the  values  of  y  for  which  the  null 
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Neighborhood  function: 

(a)  Propose  new  distribution 
parameters  using  a 
neighborhood  step 

(b)  Verify  feasibility  condition 
(3)  (i.e.,  parameters 
within  the  contour  plot) 


Calculation  of  the  new 
system’s  temperature 

_  SSEaccepted  Esbefore 


Preprocessing: 

(a)  Obtain  the  distribution  and  distribution  parameters  for  each  component 
considering  censored  data 

(b)  Obtain  the  historical  percentage  of  failures  per  component 


Termination  criteria 
determined  by 

-  Maximum  trials  ( trails ) 

-  Accepted  trials  ( trails ) 

-  Markov  chains  without 
improvement  ( chains ) 


Monte  Carlo  simulation 

(a)  Generate  failure  times 

(b)  Compute  the  system’s  time  to  failure  (TTF) 

(c)  Compute  the  simulated  percentages  of  failures 
per  component 


Single -objective  function: 
Compute  the  sum  of  square  errors  (SSE)  (2) 


SSElast  accepted 


Metropolis  criterion: 

new  parameters  are  accepted  if 

L7(0  1)  <  g_SSEnew+SSEbefore/Ts 


New  state  (set  of 
parameters) 
is  accepted 


Figure  3:  Schematic  diagram  of  the  SA-based  optimization  procedure. 


hypothesis,  y  =  y0,  is  not  rejected  at  some  stated  level  of 
significance  (1  -  a).  In  other  words,  it  consists  of  all  values 
of  y0  satisfying 

InL*  (y0)<\nL*  ~^x2i,a-  (9) 

The  LR  confidence  intervals  can  be  graphically  identified 
with  the  use  of  the  likelihood  contours,  which  consist  of 
all  the  values  of  y  and  /3  for  which  In  L(y,  /3)  is  a  constant. 
Specifically,  the  contours  constructed  are  solutions  to 

lni(rT)  =  lni*  -  \xha-  d°) 

Solutions  to  In  L*(y0)  =  In  L*  -  (1/2)^21  a  will  lie  on  these 
contours.  For  the  Weibull,  lognormal,  and  normal  distribu¬ 
tions,  the  confidence  intervals  can  be  graphically  estimated  by 
drawing  lines  that  are  both  perpendicular  to  the  coordinate 
axis  of  interest  and  tangent  to  the  likelihood  contour  [18] .  For 
one-parameter  distributions  (e.g.,  exponential),  no  contour 
plot  is  generated  and  the  confidence  bound  on  the  parameter 
will  be  considered.  The  summary  of  steps  to  obtain  the  LR 
based  bounds  on  parameters  follows. 

Step  1.  Load  the  experimental  data  and  define  separate  sets 
of  information  containing  the  failure  and  the  censored  data 
points  from  the  input  data  set. 

Step  2.  Define  the  custom  probability  density  function  (PDF) 
and  cumulative  density  function  (CDF)  equations  based  on 
the  distribution  of  the  component. 


Step  3.  Use  the  “mlecustom()”  function  in  MATLAB  to 
estimate  the  distributions  parameters  values. 

Step  4.  Compute  the  output  value  of  the  likelihood  function 
(i.e.,  the  modified  likelihood  function  for  a  component  with 
censored  data  points)  for  the  parameter  estimates  obtained  in 
Step  3  using  (5). 

Step  5.  Obtain  the  graphical  estimate  of  confidence  intervals 
of  the  parameters  satisfying  (10)  from  the  likelihood  contour 
plot. 

Step  6.  Obtain  the  chi  square  statistic  value  and  substitute 
that  value  in  (11)  [19]  for  the  specified  confidence  level: 

<n) 

where  L(6)  is  the  likelihood  function  for  the  unknown 
parameter  vector  0,  L(6)  is  the  likelihood  function  calculated 
at  the  estimated  vector  0,  and  ^ a.k  is  the  chi- squared  statistic 
with  probability  1  -  a  and  k  degrees  of  freedom,  where  k  is 
the  number  of  quantities  jointly  estimated. 

Step  7.  Using  the  graphical  estimates  of  confidence  intervals 
on  the  parameters  as  initial  guess,  two  nonlinear  optimization 
problems  are  solved  to  obtain  accurate  LR  limits.  Since  there 
are  two  unknowns  (e.g.,  Weibull s  shape  and  scale  parame¬ 
ters)  and  only  one  equation,  use  an  iterative  method  to  obtain 
the  values  of  the  parameters  (i.e.,  for  given  values  of  scale, 
obtain  the  maximum  and  minimum  value  of  shape  parameter 
that  satisfy  (11))  and  vice  versa.  The  “fsolveQ”  function  in 
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Figure  4:  Schematic  flowchart  of  the  proposed  method. 


MATLAB  was  used  to  compute  solution  iteratively.  The 
default  algorithm  for  “fsolve()”  is  the  “Trust-Region  Dogleg,” 
which  was  developed  by  Powell  [20,  21]. 

3.3.  Summary  of  Proposed  Solution  Procedure.  The  following 
steps  summarize  the  proposed  methodology. 

Step  1.  The  initial  data,  consisting  of  failure  times  and 
operating  times  of  unfailed  units,  is  provided.  The  input  data 
is  used  to  determine  the  distribution  behavior  of  the  com¬ 
ponents  failure  rate  and  to  determine  the  initial  distribution 
parameters  values  for  each  component  using  Weibull++. 

Step  2.  The  components  distributions  and  their  parameters 
are  used  to  simulate  the  expected  failure  time  ( t ).  For  each 
component,  a  random  number,  that  is,  F(t)  in  the  range  0  < 
F(t)  <  1,  is  generated  using  uniform  distribution  (7(0, 1). 
A  failure  time  is  generated  corresponding  to  F(t)  from  the 
inverse  cumulative  distribution  function  (CDF)  equation  for 
each  component.  Finally,  the  systems  time  to  failure  (TTF)  is 
computed  using  the  systems  equation  (1). 

Step  3.  The  simulated  percentages  of  failure  are  obtained 
from  the  Monte  Carlo  simulation  for  all  the  components 
in  the  system.  The  discrepancies  (e.g.,  SSE)  between  the 
simulated  and  historical  percentages  of  failures  are  computed. 

Step  4.  The  simulated  annealing  algorithm  minimizes  the 
sum  of  square  errors  (SSE)  subject  to  constraints  (e.g.,  the 
confidence  bounds  on  distribution  parameters). 

A  schematic  diagram  of  the  proposed  method  is  shown  in 
Figure  4. 


Table  2:  Distribution  fittings. 


Distribution 

Parameters 

Component  number 

Weibull  (2P) 

Shape  (jS),  scale  (a) 

1,  5,  6,  7,  8,  9, 10, 11, 12, 13, 
14, 15, 16, 19,  20,  23,  24, 

25,  26,  27,  28,  29,  30,  31, 

32,  33,  34,  35,  36,  37,  38, 

39,  40,  41,  42,  43,  45,  46, 
47,  48,  49,  50,  51,  52,  53, 

54 

Log  normal 

log (A),  log-a  (0 

2,3,17,18,22,44 

Normal 

Mean  (p),  standard 
deviation  ( o ) 

4 

Exponential 

Mean  (p) 

21 

4.  Case  Study 

4.1.  Distribution  and  Parameter  Fittings.  The  turbine  engine 
model  consists  of  54  components/failure  modes.  Table  2 
summarizes  the  parameters  corresponding  to  the  distribution 
of  each  component  used  in  the  analysis.  To  select  the 
appropriate  distribution  and  parameter  values,  a  weighted 
goodness  of  fit  measure  using  Weibull++  software  along  with 
engineering  experience  was  utilized.  A  fixed  cumulative  time 
of  1850  hours  is  used  for  the  18th  component  of  the  model. 

4.2.  The  SA’s  Parameters  Tuning  through  Statistical  Design  of 
Experiments  (DOE).  A  statistical  design  of  experiments  can 
be  defined  as  a  series  of  tests  in  which  purposeful  changes 
are  made  to  the  input  variables  of  a  process  or  system  so 
that  changes  in  the  output  response  can  be  observed  and 
quantified  [22]. 
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Table  3:  Factors  and  levels  during  the  DOE  development. 


Factor  name 

1 

Levels 

2 

3 

Trials 

50 

100 

200 

Chains 

1 

2 

4 

Alpha 

0.9 

0.95 

0.98 

Main  effects  plot  for  SSE 
Means 


Figure  5:  Main  effects  plot  for  SSE.  The  gray  background  represents 
terms  that  were  not  significant  in  the  ANOVA. 


In  order  to  tune  the  parameters  of  the  optimization  algo¬ 
rithm,  a  design  of  experiments  is  performed  to  understand 
the  cause  and  effect  relationship  among  the  SAs  parameters 
and  the  output  response. 

The  factors  that  were  selected  during  the  DOE  develop¬ 
ment  are  (1)  the  number  of  accepted  trials  before  the  systems 
temperature  is  decreased  ( Trials ),  (2)  the  number  of  Markov 
chains  (Chains),  and  (3)  the  rate  of  cooling  (Alpha).  Table  3 
shows  the  selected  factors  and  their  levels. 

The  responses  are  (1)  the  sum  of  squared  errors  (SSE),  (2) 
the  mean  of  squared  errors  (MSE),  and  (3)  the  average  CPU 
time. 

5.  Results  and  Discussion 

5.1.  DOE  Results.  The  DOE  analysis  was  performed  using 
Minitab.  A  total  number  of  33  =  27  computational  experi¬ 
ments  were  performed  while  considering  a  95%  significance 
level.  The  ANOVA  model  adequacy  verification  did  not 
exhibit  issues  with  regard  to  (1)  the  residuals  normality 
behavior,  (2)  the  residuals  independency,  and  (3)  the  resid¬ 
uals  constant  variance  assumptions. 

From  the  ANOVA  results,  it  was  observed  that  the  cooling 
rate  (Alpha)  was  not  a  significant  factor  and  it  was  fixed  to 
0.98.  The  significant  factors  were  number  of  accepted  trials 
(Trials)  and  the  number  of  Markov  chains  (Chains). 

For  the  sake  of  brevity,  only  selected  plots  are  presented. 
Figure  5  shows  the  main  effects  plot  considering  SSE  as  a 
response. 

From  Figure  5,  it  can  be  observed  that  both  factors  (e.g., 
Trials  and  Chains)  have  great  impact  on  the  response  (i.e., 


Table  4:  Comparison  of  the  responses  (before  and  after  the 
optimization). 


Sum  of  squared  errors 

Mean  of  squared 

(SSE) 

errors  (MSE) 

Before  optimization 

50.6541 

0.9380 

After  optimization 

35.1834 

0.6515 

SSE).  After  a  careful  inspection  of  the  DOE  results,  the  best 
combination  of  SAs  parameter  values  was  found  to  be 

(i)  cooling  rate  (Alpha):  0.98, 

(ii)  number  of  accepted  trials  (Trials):  200, 

(iii)  number  of  Markov  chains  (Chains):  4. 

The  initial  systems’  temperature  was  fixed  to  10.  In  order 
to  get  this  value,  the  initial  temperature  is  initialized  at  a  very 
low  value  and  an  experimental  run  is  performed  where  the 
temperature  is  raised  by  a  ratio  until  a  minimum  acceptance 
ratio  defined  by  the  practitioner  (0.8)  is  met;  the  acceptance 
ratio  is  computed  as  the  number  of  accepted  trials  (Trials) 
divided  by  the  trials. 

A  sensitivity  analysis  was  conducted  to  determine  the 
best  value  for  the  Monte  Carlo  sample  (MC  is  not  an  SA 
algorithmic  parameter  but  it  is  a  parameter  in  the  overall 
solution  procedure).  It  was  observed  that  the  SSE  becomes 
stable  for  Monte  Carlo  sample  size  of 10,000  and  larger  values; 
that  is,  the  SSE  is  not  sensible  to  changes  in  the  Monte  Carlo 
sample  size  after  a  value  of 10,000  samples.  On  the  other  hand, 
if  the  sample  size  for  the  Monte  Carlo  loop  is  below  10,000, 
the  SSE  considerably  increases  for  this  particular  problem. 
Hence,  the  selected  number  of  Monte  Carlo  samples  was  fixed 
at  10,000. 

5.2.  Test  Case  Results.  The  evaluation  of  the  test  case  was 
performed  based  on  the  outcomes  of  the  previously  presented 
DOE  analysis.  After  conducting  a  computational  evaluation 
of  the  test  case,  the  value  of  the  SSE  was  found  to  be 
35.1834  having  an  average  CPU  time  of  55,239  seconds. 
Table  4  compares  the  values  of  the  SSE  before  and  after  the 
optimization. 

The  SA-based  procedure  was  executed  five  times  con¬ 
sidering  the  best  combination  of  algorithmic  parameters. 
Figure  6  shows  the  average  performance  of  the  solution 
procedure. 

Table  5  shows  the  discrepancies  between  the  simulated 
and  historical  percentages  of  failures  for  the  components 
before  and  after  the  optimization. 

Figure  7  shows  the  histogram  of  percentage  of  failures  for 
the  historical  data,  before  and  after  the  optimization. 

In  order  to  observe  the  shifting  of  the  distributions 
parameters  at  the  different  evaluations  of  the  DOE,  contour 
plots  were  constructed  for  each  component.  Figure  8  shows 
the  shifting  of  the  distributions  parameters  at  the  different 
evaluations  of  the  DOE  for  the  45th  component  and  Figure  9 
shows  both  the  original  (fitted)  and  best  found  combination 
of  distribution  parameters  values  for  45th  component. 
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Table  5:  Comparison  of  the  responses  (before  and  after  the  optimization). 


Component  number 

Historical 
percentage 
of  failures 

Simulated 
percentage  of 
failure 

Before 

Discrepancy  [%] 

Squared  error 

Simulated 
percentage  of 
failure 

After 

Discrepancy  [%] 

Squared  error 

1 

0.86 

0.5211 

-0.3389 

0.1149 

0.46 

-0.4 

0.1600 

2 

1.82 

1.5673 

-0.2527 

0.0639 

1.66 

-0.16 

0.0256 

3 

3.53 

3.6 

0.07 

0.0049 

3.76 

0.23 

0.0529 

4 

0.88 

0.5947 

-0.2853 

0.0814 

0.63 

-0.25 

0.0625 

5 

2.4 

2.5271 

0.1271 

0.0162 

2.37 

-0.03 

0.0009 

6 

1.31 

1.2096 

-0.1004 

0.0101 

1.22 

-0.09 

0.0081 

7 

0.76 

0.6553 

-0.1047 

0.0110 

0.78 

0.02 

0.0004 

8 

0.4 

3.5392 

3.1392 

9.8546 

3.04 

2.64 

6.9696 

9 

1.04 

0.1024 

-0.9376 

0.8791 

0.09 

-0.95 

0.9025 

10 

5.35 

5.0636 

-0.2864 

0.0820 

5.09 

-0.26 

0.0676 

11 

4.56 

4.21 

-0.35 

0.1225 

4.26 

-0.3 

0.0900 

12 

8.32 

7.9168 

-0.4032 

0.1626 

7.87 

-0.45 

0.2025 

13 

2.5 

2.0153 

-0.4847 

0.2349 

2.38 

-0.12 

0.0144 

14 

1.39 

1.1423 

-0.2477 

0.0614 

1.18 

-0.21 

0.0441 

15 

2.8 

2.5067 

-0.2933 

0.0860 

2.45 

-0.35 

0.1225 

16 

0.3 

0.259 

-0.041 

0.0017 

0.27 

-0.03 

0.0009 

17 

0.51 

0.499 

-0.011 

0.0001 

0.47 

-0.04 

0.0016 

18 

0.91 

1.8601 

0.9501 

0.9027 

1.81 

0.9 

0.8100 

19 

1.34 

0.9799 

-0.3601 

0.1297 

1.24 

-0.1 

0.0100 

20 

0.68 

0.4956 

-0.1844 

0.0340 

0.68 

0 

0.0000 

21 

0.45 

0.4621 

0.0121 

0.0001 

0.53 

0.08 

0.0064 

22 

0.51 

0.4554 

-0.0546 

0.0030 

0.4 

-0.11 

0.0121 

23 

0.25 

0.1628 

-0.0872 

0.0076 

0.17 

-0.08 

0.0064 

24 

11.07 

6.3014 

-4.7686 

22.7395 

7.03 

-4.04 

16.3216 

25 

4.03 

4.1075 

0.0775 

0.0060 

4.26 

0.23 

0.0529 

26 

0.55 

0.0571 

-0.4929 

0.2430 

0.05 

-0.5 

0.2500 

27 

0.7 

0.6549 

-0.0451 

0.0020 

0.72 

0.02 

0.0004 

28 

5.8 

4.4421 

-1.3579 

1.8439 

4.31 

-1.49 

2.2201 

29 

0.38 

0.2823 

-0.0977 

0.0095 

0.22 

-0.16 

0.0256 

30 

3.32 

2.1621 

-1.1579 

1.3407 

2.35 

-0.97 

0.9409 

31 

0.5 

0.4702 

-0.0298 

0.0009 

0.43 

-0.07 

0.0049 

32 

0.99 

0.7867 

-0.2033 

0.0413 

0.91 

-0.08 

0.0064 

33 

0.45 

0.3763 

-0.0737 

0.0054 

0.33 

-0.12 

0.0144 

34 

0.4 

0.3421 

-0.0579 

0.0034 

0.38 

-0.02 

0.0004 

35 

0.63 

0.822 

0.192 

0.0369 

0.81 

0.18 

0.0324 

36 

0.85 

1.0452 

0.1952 

0.0381 

0.98 

0.13 

0.0169 

37 

0.4 

0.561 

0.161 

0.0259 

0.47 

0.07 

0.0049 

38 

0.38 

0.5465 

0.1665 

0.0277 

0.49 

0.11 

0.0121 

39 

0.23 

0.3274 

0.0974 

0.0095 

0.36 

0.13 

0.0169 

40 

6.02 

6.3466 

0.3266 

0.1067 

6.39 

0.37 

0.1369 

41 

2.98 

4.0667 

1.0867 

1.1809 

3.88 

0.9 

0.8100 

42 

0.61 

0.9346 

0.3246 

0.1054 

0.93 

0.32 

0.1024 

43 

0.6 

0.8256 

0.2256 

0.0509 

0.91 

0.31 

0.0961 

44 

0.76 

0.2627 

-0.4973 

0.2473 

0.27 

-0.49 

0.2401 
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Table  5:  Continued. 


Component  number 

Historical 
percentage 
of  failures 

Simulated 
percentage  of 
failure 

Before 

Discrepancy  [%] 

Squared  error 

Simulated 
percentage  of 
failure 

After 

Discrepancy  [%] 

Squared  error 

45 

6.28 

9.0659 

2.7859 

7.7612 

7.84 

1.56 

2.4336 

46 

1.23 

1.8439 

0.6139 

0.3769 

1.55 

0.32 

0.1024 

47 

0.85 

1.1438 

0.2938 

0.0863 

1.1 

0.25 

0.0625 

48 

1.66 

2.1355 

0.4755 

0.2261 

2.37 

0.71 

0.5041 

49 

1.29 

1.8673 

0.5773 

0.3333 

1.76 

0.47 

0.2209 

50 

0.71 

0.9715 

0.2615 

0.0684 

1.01 

0.3 

0.0900 

51 

2.11 

2.9768 

0.8668 

0.7513 

2.92 

0.81 

0.6561 

52 

0.38 

0.558 

0.178 

0.0317 

0.68 

0.3 

0.0900 

53 

0.5 

0.7342 

0.2342 

0.0548 

0.82 

0.32 

0.1024 

54 

0.45 

0.6368 

0.1868 

0.0349 

0.66 

0.21 

0.0441 

Total 

50.6541 

Total 

35.1834 

-  Average  +  Std  Dev  -  Average  -  Std  Dev 

- Average 


Figure  6:  The  average  performance  of  the  SA-based  procedure. 


■  Historical  percentage  of  failures 

■  Simulated  percentage  of  failures  (before) 

■  Simulated  percentage  of  failures  (after) 

Figure  7:  Histogram  of  percentages  of  failures  before  and  after  optimization. 
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-*■  Distributions  parameters  for  different  DOE  evaluations 
^  Distributions  parameter  using  the  tuned  SAs  parameters 


Figure  8:  Distribution  parameters  for  the  45th  component  for 
different  evaluation  cases. 


^  Original  (fitted)  distributions  parameters  values 
^  Best  found  distributions  parameters  values 


Figure  9:  Original  distributions  parameters  values  and  best  found 
combination  of  distributions  parameters  values  for  the  45th  compo¬ 
nent. 

6.  Conclusions  and  Future  Work 

This  paper  presented  a  simulated  annealing-based  optimiza¬ 
tion  method  to  minimize  the  discrepancy  of  historical  and 
simulated  percentages  of  failures  in  a  turbine  engine  model. 
A  DOE  was  performed  for  the  tuning  of  the  algorithmic 


parameters.  The  results  showed  a  30%  reduction  of  the 
SSE  (e.g.,  from  50.6541  to  35.1834)  for  the  engine  model. 
The  average  CPU  time  was  approximately  15  hours  mainly 
due  to  the  calculations  involved  in  the  likelihood  function. 
Alternative  neighborhood  and  feasibility  functions  can  be 
investigated  by  studying  the  trends  in  the  shifting  parameters 
values  per  component.  For  instance,  it  was  observed  that 
the  distributions  parameters  shifted  around  the  edge  of  the 
contour  plot  for  some  components. 

The  proposed  simulation-based  optimization  method  can 
serve  as  a  decision-making  tool  for  maintenance,  repair, 
and  overhaul  companies  and  will  potentially  reduce  the  cost 
of  labor  associated  to  finding  the  appropriate  value  of  the 
distributions  parameters  for  each  component/failure  mode  in 
the  model  and  increase  the  accuracy  in  the  prediction  of  the 
mean  time  to  failures  (MTTF). 

Future  research  lines  involve  parallelization  of  the  algo¬ 
rithm  to  solve  larger  models  (e.g.,  thousands  of  components) 
and  comparing  the  performance  of  the  simulated  annealing 
with  other  metaheuristics  such  as  Evolutionary  Algorithms, 
Tabu  Search,  and  Particle  Swarm  Optimization,  among  oth¬ 
ers. 

Appendix 

Error  Function.  The  error  function  erf(X)  is  twice  the  integral 
of  the  Gaussian  distribution  with  0  mean  and  variance  of  1/2. 
The  equation  of  error  function  is 

erf  (x)  =  [  e~x  dt.  (A.l) 

V71  Jo 

The  complementary  error  function  or  erfc(X)  is  the 
complement  of  the  error  function;  that  is,  erfc(X)  =  1  - 
erf(X) 

2  f  ^  2 

erfc  (x)  =  ——  e~x  dt.  (A.2) 

V n  Jx 

Nomenclature 

MTTF:  Mean  time  to  failure  of  the  system 
STTF:  Standard  deviation  of  time  to  failure  of  the 
system 

TTF:  Time  to  failure  of  the  system 

g(TK ):  System  TTF 

Tk:  Vector  of  K  variables 

t:  Failure  time 

F(t):  Cumulative  distribution  function  of 

failure  time  (f ) 

n:  Number  of  components  in  the  engine 

model 

a:  Scale  parameter  of  a  Weibull  distribution 

/?:  Shape  parameter  of  a  Weibull  distribution 

A:  Mean  of  log(x)  of  a  lognormal  distribution 

C:  Standard  deviation  of  log(x)  of  a 

lognormal  distribution 
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[i\  Mean  of  a  normal  distribution  or  an 
exponential  distribution 
a:  Standard  deviation  of  a  normal 
distribution. 
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